//construct provider propensity without controls. Include hospital-time FE only. 


use "$data_out/naive.dta", clear

gen propensity = .

gen month = mofd(encdate)

gen dow = dow(encdate)
egen hym  = group(month clinic)
egen hd = group(clinic dow)
gen temp_dx3 = substr(dx1,1,3)
egen diag = group(temp_dx3)



*drop missing and singleton observations that won't be included in regressions
keep if !mi(age_bin) & !mi(mos) & !mi(longevity) & !mi(married) & !mi(junior) & !mi(college) & !mi(race_white) & !mi(female) & !mi(diag) & !mi(afqt_p) 
gen include = 1
foreach i in diag mos hym hd {
bys `i': replace include = 0 if _N == 1
}
keep if include == 1









reghdfe ed_prescription , absorb(hym hd) res(residual)
egen mean_residual = mean(residual), by(provID year)

bys provID year: gen N = _N
cap drop propensity
gen propensity = (mean_res * N - residual)/(N-1)
drop N
sort provID year


preserve 
bys provID year: keep if _n ==1

bys year (mean_residual): gen n = _n 
bys year (mean_residual): gen N = _N 
gen above_median = n > N/2
gen quartile = 1 if n <= N/4
replace quartile = 2 if n <= N/2 & n > N/4
replace quartile = 3 if n > N/2 & n <= 3*N/4
replace quartile = 4 if n > 3*N/4
keep provID year above_median quartile

tempfile a
save `a'
restore

merge m:1 provID year using `a', nogen







save "$data_out/provider_propensity_naive1", replace



//add in diagnosis code


 use "$data_out/naive.dta", clear

gen propensity = .

gen month = mofd(encdate)

gen dow = dow(encdate)
egen hym  = group(month clinic)
egen hd = group(clinic dow)
gen temp_dx3 = substr(dx1,1,3)
egen diag = group(temp_dx3)



*drop missing and singleton observations that won't be included in regressions
keep if !mi(age_bin) & !mi(mos) & !mi(longevity) & !mi(married) & !mi(junior) & !mi(college) & !mi(race_white) & !mi(female) & !mi(diag) & !mi(afqt_p) 
gen include = 1
foreach i in diag mos hym hd {
bys `i': replace include = 0 if _N == 1
}
keep if include == 1





use "$data_out/provider_propensity_naive2", clear


*merge indicators for depresson
merge 1:m pid_pde encdate using "$data_out/depression_scripts",nogen keep(1 3)
merge 1:m pid_pde encdate using "$data_out/APFT_pre", nogen keep(1 3)
merge 1:1 pid_pde encdate using "$data_out/deployments_pre", nogen keep(1 3)
merge 1:1 pid_pde encdate using "$data_out/Clearances_pre", nogen keep(1 3)
merge 1:1 pid_pde encdate using "$data_out/Drug_Tests_Pre", nogen keep(1 3)
merge 1:1 pid_pde encdate using "$data_out/Flags_pre", nogen keep(1 3)
merge 1:1 pid_pde encdate using "$data_out/Promotions_Pre", nogen  keep(1 3)
merge 1:1 pid_pde encdate using "$data_out/PULHES_pre", nogen keep(1 3)
merge 1:1 pid_pde encdate using "$data_out/depression_pre", nogen keep(1 3)


rename (depression_diagnosis subsample) (depression_pre depression_post)
recode depression_post .=0
recode depression_pre .=0

label var propensity "Physician Propensity to Prescribe Opioids"
label var age "Age"
label var age_bin "Age Bin"
	label define age_labels 10 "Under 20" 20 "20-39" 30 "30-39" 40 "40-49" 50 "50-59" 60 "60-69" 70 "70+"
	label values age_bin age_labels
label var race_white "Race - White"
label var college "Education - College Degree"
label var female "Female"
label var junior_enlisted "Rank - Junior"
label var married "Marital Status - Married"
label var longevity "Military Tenure (Yrs)"
label var depression_post "Depression or Anxiety"
label var afqt_pct "Armed Forces Qualificaiton Test Percentile"



			cap postclose balance
			postfile balance str16 Variable Mean Coefficient Std_Error p_value   using "$plotdir/balance_continuous_naive2.dta", replace
			
			local dvar "z_score_pre deploy_pre denied_pre revoke_pre physical_pre_Z psych_pre_Z any_flag_pre substance_flag_pre  crime_flag_pre adverse_flag_pre promotion_pre demotion_pre depression_pre depression_post  age female race_white junior_enlisted college married longevity afqt_p" 

foreach var of varlist `dvar' { 
			eststo balance: reg `var'  propensity ,  cluster(provID)
			estadd ysumm, replace
		
			post balance ("`var'") (e(ymean)) (_b[propensity]) (_se[propensity]) (2*normal(-abs(_b[propensity]/_se[propensity])))
			
			}

	postclose balance
	
	*create indicator to drop singletons for binary instrument in summary stats

preserve 
keep if inlist(quartile,1,4)
cap drop include
foreach i in diag mos hym hd {
bys `i': replace include = 0 if _N == 1
}
keep if include == 1

			local dvar "z_score_pre deploy_pre denied_pre revoke_pre physical_pre_Z psych_pre_Z any_flag_pre substance_flag_pre  crime_flag_pre adverse_flag_pre promotion_pre demotion_pre depression_pre depression_post  age female race_white junior_enlisted college married longevity afqt_p" 

	
				cap postclose balance
			postfile balance str16 Variable  Coefficient Std_Error p_value   using "$plotdir/balance_binary_naive2.dta", replace

foreach var of varlist `dvar' { 
	eststo balance: reg `var'  above  if inlist(quartile,1,4), cluster(provID)
			
		
			post balance ("`var'")  (_b[above]) (_se[above]) (2*normal(-abs(_b[above]/_se[above])))
			
			}

	postclose balance
	restore


			cap postclose balance
			postfile balance str16 Variable Mean Coefficient Std_Error p_value   using "$plotdir/balance_continuousFE_naive`i'.dta", replace

foreach var of varlist  age female race_white junior_enlisted college married longevity depression afqt_p { 
			eststo balance: reghdfe `var'  propensity , absorb(hym hd) cluster(provID)
			estadd ysumm, replace
		
			post balance ("`var'") (e(ymean)) (_b[propensity]) (_se[propensity]) (2*normal(-abs(_b[propensity]/_se[propensity])))
			
			}

	postclose balance
	
				cap postclose balance
			postfile balance str16 Variable  Coefficient Std_Error p_value   using "$plotdir/balance_binaryFE_naive`i'.dta", replace

foreach var of varlist  age female race_white junior_enlisted college married longevity depression  afqt_p { 
	eststo balance: reghdfe `var'  above  if inlist(quartile,1,4), absorb(hym hd) cluster(provID)
			
		
			post balance ("`var'")  (_b[above]) (_se[above]) (2*normal(-abs(_b[above]/_se[above])))
			
			}

	postclose balance





//table
use "$plotdir/balance_continuous_naive2.dta", clear
rename ( Coefficient Std_Error p_value) (Continuous_Coef Continuous_SE Continuous_P)
merge 1:1 Variable using "$plotdir/balance_binary_naive2.dta", nogen 
rename ( Coefficient Std_Error p_value) (Binary_Coef Binary_SE Binary_P)
export excel using "\\facaa7pde01p\PDE\Projects\ACADEMIC\ACADEMIC-UOA\Draft - July 2021\Revisions - AEJ Applied\Tables and Figures\balance_naive.xls", sheetreplace firstrow(variables)






eststo balance_script: reg ed_prescription age race_white female junior_enlisted college married longevity depression afqt_p,  cluster(provID)

eststo balance_propensity: reg propensity age race_white female junior_enlisted college married longevity depression afqt_p, cluster(provID)


#delimit ;
	//word
esttab balance_script balance_propensity using "EZ_Table.rtf" , 
	replace label se star(* 0.10 ** 0.05 *** 0.01)
	b(4) se(4)
	stats(N F, 
	fmt(%-9.0fc %-9.2fc)
	label("Observations" "F-Statistic"))
	mtitles("ED Prescription" "Residualized Instrument")
	noconstant nobaselevels
	title(Balance Table for Continuous Instrument)
	compress;
	

#delimit cr

//eichmeyer and zhang table

forval i = 1/2 {
use "$data_out/provider_propensity_naive`i'", clear


*merge indicators for depresson
merge 1:m pid_pde encdate using "$data_out/depression_pre",nogen keep(1 3)


recode depression .=0

label var propensity "Physician Propensity to Prescribe Opioids"
label var age "Age"
label var age_bin "Age Bin"
	label define age_labels 10 "Under 20" 20 "20-39" 30 "30-39" 40 "40-49" 50 "50-59" 60 "60-69" 70 "70+"
	label values age_bin age_labels
label var race_white "Race - White"
label var college "Education - College Degree"
label var female "Female"
label var junior_enlisted "Rank - Junior"
label var married "Marital Status - Married"
label var longevity "Military Tenure (Yrs)"
label var depression "Depression or Anxiety"
label var afqt_pct "Armed Forces Qualificaiton Test Percentile"


eststo balance_script: reg ed_prescription age race_white female junior_enlisted college married longevity depression afqt_p,  cluster(provID)

eststo balance_propensity: reg propensity age race_white female junior_enlisted college married longevity depression afqt_p, cluster(provID)


#delimit ;
	//word
esttab balance_script balance_propensity using "EZ_Table_naive`i'.rtf" , 
	replace label se star(* 0.10 ** 0.05 *** 0.01)
	b(4) se(4)
	stats(N F, 
	fmt(%-9.0fc %-9.2fc)
	label("Observations" "F-Statistic"))
	mtitles("ED Prescription" "Residualized Instrument")
	noconstant nobaselevels
	title(Balance Table for Continuous Instrument)
	compress;
	

#delimit cr

}

